
#################################################################################################

#  Description: Replicating robustness checks (Appendix B) from Fleming, Thomas G. (Forthcoming) 
#               'Why Change a Winning Team? Explaining Post-Election Cabinet Reshuffles in Four
#                Westminster Democracies', Political Studies.
#  Author:      Thomas G. Fleming
#  Date:        10/09/2021
#  Note:        Requires file "reshuffles_data.csv" to be saved in same folder

#################################################################################################

rm(list=ls())

library(stargazer)  # for stargazer()
library(MASS)       # for glm.nb()
library(DescTools)  # for PseudoR2

library(rstudioapi) # sets location as wd
current_path <- getActiveDocumentContext()$path
setwd(dirname(current_path))

####################################
# Loads data and formats variables #
####################################

  data <- read.csv("./reshuffles_data.csv")
  
  data$election.date <- as.Date(data$election.date, format="%Y-%m-%d") # formats as date
  data$cabinet.date <- as.Date(data$cabinet.date, format="%Y-%m-%d")   # formats as date
  
  data$selection.rule <- relevel(data$selection.rule, ref = "pm")      # sets reference category
  
  # generating logged version of seat share
  
  data$ln.seatshare <- log(data$seatshare)
  
################################################################################################
# Creates dependent variable - number of pre-election ministers moved within or out of cabinet #
################################################################################################

  data$cabinet.changes <- data$moved + data$left
    
  # generating broader measure of dependent variable
  
  data$cabinet.changes.broad <- data$moved + data$left + data$entered
      
#####################################################################################
# Table B1: descriptive statistics for additional variables not reported in Table 1 #
#####################################################################################
  
  varlist <- data[c("cabinet.changes.broad", "voteshare", "voteshare.change", "ln.seatshare")]
  
  lapply(varlist, mean, na.rm = TRUE)
  lapply(varlist, sd, na.rm = TRUE)
  lapply(varlist, min, na.rm = TRUE)
  lapply(varlist, max, na.rm = TRUE)
  
###########################################################################################
# Table B2: Repeats original Table 2 with alternative calculation of dependent variable  #
###########################################################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.bivariate.widerDV <- glm.nb(cabinet.changes.broad ~
                                          seatshare +
                                          preelection.total,
                                        data = data)
  seatshare.bivariate.widerDV.R2 <- PseudoR2(seatshare.bivariate.widerDV, which = "McFadden")
  
  # adding various control variables
  
  seatshare.controls.widerDV <- glm.nb(cabinet.changes.broad ~
                                         seatshare +
                                         preelection.total +
                                         unavailable +
                                         median.post.tenure +
                                         median.cabinet.tenure +
                                         new.PM +
                                         selection.rule,
                                       data = data)
  seatshare.controls.widerDV.R2 <- PseudoR2(seatshare.controls.widerDV, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.countryFE.widerDV <- glm.nb(cabinet.changes.broad ~
                                          seatshare +
                                          preelection.total +
                                          unavailable +
                                          median.post.tenure +
                                          median.cabinet.tenure +
                                          new.PM +
                                          selection.rule +
                                          as.factor(country),
                                        data = data)
  seatshare.countryFE.widerDV.R2 <- PseudoR2(seatshare.countryFE.widerDV, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.partyFE.widerDV <- glm.nb(cabinet.changes.broad ~
                                        seatshare +
                                        preelection.total +
                                        unavailable +
                                        median.post.tenure +
                                        median.cabinet.tenure +
                                        new.PM +
                                        selection.rule +
                                        as.factor(party),
                                      data = data)
  seatshare.partyFE.widerDV.R2 <- PseudoR2(seatshare.partyFE.widerDV, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("seatshare", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "selection.rulecaucus", 
           "Constant")
  
  stargazer(seatshare.bivariate.widerDV, 
            seatshare.controls.widerDV, 
            seatshare.countryFE.widerDV, 
            seatshare.partyFE.widerDV,
            type = "html",
            title = "Table B2. Negative binomial regression of post-election cabinet changes",
            out = "./tableB2.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Seat share", "Total ministers", "Unavailable ministers", "Median tenure (current post)", "Median tenure (any post)",
                                 "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved or added",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B1)", "(B2)", "(B3)", "(B4)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.bivariate.widerDV), nobs(seatshare.controls.widerDV), 
                               nobs(seatshare.countryFE.widerDV), nobs(seatshare.partyFE.widerDV)),
                             c("Pseudo R2", round(seatshare.bivariate.widerDV.R2, digits = 3),
                               round(seatshare.controls.widerDV.R2, digits = 3),
                               round(seatshare.countryFE.widerDV.R2, digits = 3),
                               round(seatshare.partyFE.widerDV.R2, digits = 3))),
            font.size = "small")

##########################################################################################
# Table B3: Repeats original Table 3 with alternative calculation of dependent variable  #
##########################################################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.change.bivariate.widerDV <- glm.nb(cabinet.changes.broad ~
                                                 seatshare.change +
                                                 preelection.total,
                                               data = data)
  seatshare.change.bivariate.widerDV.R2 <- PseudoR2(seatshare.change.bivariate.widerDV, which = "McFadden")
  
  # adding various control variables
  
  seatshare.change.controls.widerDV <- glm.nb(cabinet.changes.broad ~
                                                seatshare.change +
                                                preelection.total +
                                                unavailable +
                                                median.post.tenure +
                                                median.cabinet.tenure +
                                                new.PM +
                                                selection.rule,
                                              data = data)
  seatshare.change.controls.widerDV.R2 <- PseudoR2(seatshare.change.controls.widerDV, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.change.countryFE.widerDV <- glm.nb(cabinet.changes.broad ~
                                                 seatshare.change +
                                                 preelection.total +
                                                 unavailable +
                                                 median.post.tenure +
                                                 median.cabinet.tenure +
                                                 new.PM +
                                                 selection.rule +
                                                 as.factor(country),
                                               data = data)
  seatshare.change.countryFE.widerDV.R2 <- PseudoR2(seatshare.change.countryFE.widerDV, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.change.partyFE.widerDV <- glm.nb(cabinet.changes.broad ~
                                               seatshare.change +
                                               preelection.total +
                                               unavailable +
                                               median.post.tenure +
                                               median.cabinet.tenure +
                                               new.PM +
                                               selection.rule +
                                               as.factor(party),
                                             data = data)
  seatshare.change.partyFE.widerDV.R2 <- PseudoR2(seatshare.change.partyFE.widerDV, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("seatshare.change", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "selection.rulecaucus", 
           "Constant")
  
  stargazer(seatshare.change.bivariate.widerDV, 
            seatshare.change.controls.widerDV, 
            seatshare.change.countryFE.widerDV, 
            seatshare.change.partyFE.widerDV,
            type = "html",
            title = "Table B3. Negative binomial regression of post-election cabinet changes",
            out = "./tableB3.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Seat share change", "Total ministers", "Unavailable ministers", "Median tenure (current post)", 
                                 "Median tenure (any post)", "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved or added",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B5)", "(B6)", "(B7)", "(B8)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.change.bivariate.widerDV), nobs(seatshare.change.controls.widerDV), 
                               nobs(seatshare.change.countryFE.widerDV), nobs(seatshare.change.partyFE.widerDV)),
                             c("Pseudo R2", round(seatshare.change.bivariate.widerDV.R2, digits = 3),
                               round(seatshare.change.controls.widerDV.R2, digits = 3),
                               round(seatshare.change.countryFE.widerDV.R2, digits = 3),
                               round(seatshare.change.partyFE.widerDV.R2, digits = 3))),
            font.size = "small")
  
#######################################################################################
# Table B4: Repeats original Table 2 with vote share as the main independent variable #
#######################################################################################
  
  # bivariate (controlling just for total ministers)
  
  voteshare.bivariate <- glm.nb(cabinet.changes ~
                                  voteshare +
                                  preelection.total,
                                data = data)
  voteshare.bivariate.R2 <- PseudoR2(voteshare.bivariate, which = "McFadden")
  
  # adding various control variables
  
  voteshare.controls <- glm.nb(cabinet.changes ~
                                 voteshare +
                                 preelection.total +
                                 unavailable +
                                 median.post.tenure +
                                 median.cabinet.tenure +
                                 new.PM +
                                 selection.rule,
                               data = data)
  voteshare.controls.R2 <- PseudoR2(voteshare.controls, which = "McFadden")
  
  # adding country fixed effects
  
  voteshare.countryFE <- glm.nb(cabinet.changes ~
                                  voteshare +
                                  preelection.total +
                                  unavailable +
                                  median.post.tenure +
                                  median.cabinet.tenure +
                                  new.PM +
                                  selection.rule +
                                  as.factor(country),
                                data = data)
  voteshare.countryFE.R2 <- PseudoR2(voteshare.countryFE, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  voteshare.partyFE <- glm.nb(cabinet.changes ~
                                voteshare +
                                preelection.total +
                                unavailable +
                                median.post.tenure +
                                median.cabinet.tenure +
                                new.PM +
                                selection.rule +
                                as.factor(party),
                              data = data)
  voteshare.partyFE.R2 <- PseudoR2(voteshare.partyFE, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("voteshare", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "selection.rulecaucus", 
           "Constant")
  
  stargazer(voteshare.bivariate, 
            voteshare.controls, 
            voteshare.countryFE, 
            voteshare.partyFE,
            type = "html",
            title = "Table B4. Negative binomial regression of post-election cabinet changes",
            out = "./tableB4.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Vote share", "Total ministers", "Unavailable ministers", "Median tenure (current post)", "Median tenure (any post)",
                                 "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B9)", "(B10)", "(B11)", "(B12)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(voteshare.bivariate), nobs(voteshare.controls), 
                               nobs(voteshare.countryFE), nobs(voteshare.partyFE)),
                             c("Pseudo R2", round(voteshare.bivariate.R2, digits = 3),
                               round(voteshare.controls.R2, digits = 3),
                               round(voteshare.countryFE.R2, digits = 3),
                               round(voteshare.partyFE.R2, digits = 3))),
            font.size = "small")
  
#################################################################################################
# Table B5: Repeats original Table 3 with change in vote share as the main independent variable #
#################################################################################################
  
  # bivariate (controlling just for total ministers)
  
  voteshare.change.bivariate <- glm.nb(cabinet.changes ~
                                         voteshare.change +
                                         preelection.total,
                                       data = data)
  voteshare.change.bivariate.R2 <- PseudoR2(voteshare.change.bivariate, which = "McFadden")
  
  # adding various control variables
  
  voteshare.change.controls <- glm.nb(cabinet.changes ~
                                        voteshare.change +
                                        preelection.total +
                                        unavailable +
                                        median.post.tenure +
                                        median.cabinet.tenure +
                                        new.PM +
                                        selection.rule,
                                      data = data)
  voteshare.change.controls.R2 <- PseudoR2(voteshare.change.controls, which = "McFadden")
  
  # adding country fixed effects
  
  voteshare.change.countryFE <- glm.nb(cabinet.changes ~
                                         voteshare.change +
                                         preelection.total +
                                         unavailable +
                                         median.post.tenure +
                                         median.cabinet.tenure +
                                         new.PM +
                                         selection.rule +
                                         as.factor(country),
                                       data = data)
  voteshare.change.countryFE.R2 <- PseudoR2(voteshare.change.countryFE, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  voteshare.change.partyFE <- glm.nb(cabinet.changes ~
                                       voteshare.change +
                                       preelection.total +
                                       unavailable +
                                       median.post.tenure +
                                       median.cabinet.tenure +
                                       new.PM +
                                       selection.rule +
                                       as.factor(party),
                                     data = data)
  voteshare.change.partyFE.R2 <- PseudoR2(voteshare.change.partyFE, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("voteshare.change", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "selection.rulecaucus", 
           "Constant")
  
  stargazer(voteshare.change.bivariate, 
            voteshare.change.controls, 
            voteshare.change.countryFE, 
            voteshare.change.partyFE,
            type = "html",
            title = "Table B5. Negative binomial regression of post-election cabinet changes",
            out = "./tableB5.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Vote share change", "Total ministers", "Unavailable ministers", "Median tenure (current post)", 
                                 "Median tenure (any post)", "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B13)", "(B14)", "(B15)", "(B16)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(voteshare.change.bivariate), nobs(voteshare.change.controls), 
                               nobs(voteshare.change.countryFE), nobs(voteshare.change.partyFE)),
                             c("Pseudo R2", round(voteshare.change.bivariate.R2, digits = 3),
                               round(voteshare.change.controls.R2, digits = 3),
                               round(voteshare.change.countryFE.R2, digits = 3),
                               round(voteshare.change.partyFE.R2, digits = 3))),
            font.size = "small")
  
#############################################################
# Table B6: Repeats original Table 2 with seat share logged #
#############################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.bivariate.lnIV <- glm.nb(cabinet.changes ~
                                           ln.seatshare +
                                           preelection.total,
                                         data = data)
  seatshare.bivariate.lnIV.R2 <- PseudoR2(seatshare.bivariate.lnIV, which = "McFadden")
  
  # adding various control variables
  
  seatshare.controls.lnIV <- glm.nb(cabinet.changes ~
                                          ln.seatshare +
                                          preelection.total +
                                          unavailable +
                                          median.post.tenure +
                                          median.cabinet.tenure +
                                          new.PM +
                                          selection.rule,
                                        data = data)
  seatshare.controls.lnIV.R2 <- PseudoR2(seatshare.controls.lnIV, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.countryFE.lnIV <- glm.nb(cabinet.changes ~
                                           ln.seatshare +
                                           preelection.total +
                                           unavailable +
                                           median.post.tenure +
                                           median.cabinet.tenure +
                                           new.PM +
                                           selection.rule +
                                           as.factor(country),
                                         data = data)
  seatshare.countryFE.lnIV.R2 <- PseudoR2(seatshare.countryFE.lnIV, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.partyFE.lnIV <- glm.nb(cabinet.changes ~
                                         ln.seatshare +
                                         preelection.total +
                                         unavailable +
                                         median.post.tenure +
                                         median.cabinet.tenure +
                                         new.PM +
                                         selection.rule +
                                         as.factor(party),
                                       data = data)
  seatshare.partyFE.lnIV.R2 <- PseudoR2(seatshare.partyFE.lnIV, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("ln.seatshare", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "selection.rulecaucus", 
           "Constant")
  
  stargazer(seatshare.bivariate.lnIV, 
            seatshare.controls.lnIV, 
            seatshare.countryFE.lnIV, 
            seatshare.partyFE.lnIV,
            type = "html",
            title = "Table B6. Negative binomial regression of post-election cabinet changes",
            out = "./tableB6.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Log seat share", "Total ministers", "Unavailable ministers", "Median tenure (current post)", 
                                 "Median tenure (any post)", "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B17)", "(B18)", "(B19)", "(B20)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.bivariate.lnIV), nobs(seatshare.controls.lnIV), 
                               nobs(seatshare.countryFE.lnIV), nobs(seatshare.partyFE.lnIV)),
                             c("Pseudo R2", round(seatshare.bivariate.lnIV.R2, digits = 3),
                               round(seatshare.controls.lnIV.R2, digits = 3),
                               round(seatshare.countryFE.lnIV.R2, digits = 3),
                               round(seatshare.partyFE.lnIV.R2, digits = 3))),
            font.size = "small")
  
##############################################################################
# Table B7: Repeats original Table 2 without caucus-based cabinet selections #
##############################################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.bivariate.nocaucus <- glm.nb(cabinet.changes ~
                                           seatshare +
                                           preelection.total,
                                         data = data[data$selection.rule == "pm",])
  seatshare.bivariate.nocaucus.R2 <- PseudoR2(seatshare.bivariate.nocaucus, which = "McFadden")
  
  # adding various control variables
  
  seatshare.controls.nocaucus <- glm.nb(cabinet.changes ~
                                          seatshare +
                                          preelection.total +
                                          unavailable +
                                          median.post.tenure +
                                          median.cabinet.tenure +
                                          new.PM,
                                        data = data[data$selection.rule == "pm",])
  seatshare.controls.nocaucus.R2 <- PseudoR2(seatshare.controls.nocaucus, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.countryFE.nocaucus <- glm.nb(cabinet.changes ~
                                           seatshare +
                                           preelection.total +
                                           unavailable +
                                           median.post.tenure +
                                           median.cabinet.tenure +
                                           new.PM +
                                           as.factor(country),
                                         data = data[data$selection.rule == "pm",])
  seatshare.countryFE.nocaucus.R2 <- PseudoR2(seatshare.countryFE.nocaucus, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.partyFE.nocaucus <- glm.nb(cabinet.changes ~
                                         seatshare +
                                         preelection.total +
                                         unavailable +
                                         median.post.tenure +
                                         median.cabinet.tenure +
                                         new.PM +
                                         as.factor(party),
                                       data = data[data$selection.rule == "pm",])
  seatshare.partyFE.nocaucus.R2 <- PseudoR2(seatshare.partyFE.nocaucus, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("seatshare", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "Constant")
  
  stargazer(seatshare.bivariate.nocaucus, 
            seatshare.controls.nocaucus, 
            seatshare.countryFE.nocaucus, 
            seatshare.partyFE.nocaucus,
            type = "html",
            title = "Table B7. Negative binomial regression of post-election cabinet changes",
            out = "./tableB7.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Seat share", "Total ministers", "Unavailable ministers", "Median tenure (current post)", "Median tenure (any post)",
                                 "PM's first election", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B21)", "(B22)", "(B23)", "(B24)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.bivariate.nocaucus), nobs(seatshare.controls.nocaucus), 
                               nobs(seatshare.countryFE.nocaucus), nobs(seatshare.partyFE.nocaucus)),
                             c("Pseudo R2", round(seatshare.bivariate.nocaucus.R2, digits = 3),
                               round(seatshare.controls.nocaucus.R2, digits = 3),
                               round(seatshare.countryFE.nocaucus.R2, digits = 3),
                               round(seatshare.partyFE.nocaucus.R2, digits = 3))),
            font.size = "small")
  
##############################################################################
# Table B8: Repeats original Table 3 without caucus-based cabinet selections #
##############################################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.change.bivariate.nocaucus <- glm.nb(cabinet.changes ~
                                                  seatshare.change +
                                                  preelection.total,
                                                data = data[data$selection.rule == "pm",])
  seatshare.change.bivariate.nocaucus.R2 <- PseudoR2(seatshare.change.bivariate.nocaucus, which = "McFadden")
  
  # adding various control variables
  
  seatshare.change.controls.nocaucus <- glm.nb(cabinet.changes ~
                                                 seatshare.change +
                                                 preelection.total +
                                                 unavailable +
                                                 median.post.tenure +
                                                 median.cabinet.tenure +
                                                 new.PM,
                                               data = data[data$selection.rule == "pm",])
  seatshare.change.controls.nocaucus.R2 <- PseudoR2(seatshare.change.controls.nocaucus, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.change.countryFE.nocaucus <- glm.nb(cabinet.changes ~
                                                  seatshare.change +
                                                  unavailable +
                                                  median.post.tenure +
                                                  median.cabinet.tenure +
                                                  new.PM +
                                                  preelection.total +
                                                  as.factor(country),
                                                data = data[data$selection.rule == "pm",])
  seatshare.change.countryFE.nocaucus.R2 <- PseudoR2(seatshare.change.countryFE.nocaucus, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.change.partyFE.nocaucus <- glm.nb(cabinet.changes ~
                                                seatshare.change +
                                                preelection.total +
                                                unavailable +
                                                median.post.tenure +
                                                median.cabinet.tenure +
                                                new.PM +
                                                as.factor(party),
                                              data = data[data$selection.rule == "pm",])
  seatshare.change.partyFE.nocaucus.R2 <- PseudoR2(seatshare.change.partyFE.nocaucus, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("seatshare.change", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", "Constant")
  
  stargazer(seatshare.change.bivariate.nocaucus, 
            seatshare.change.controls.nocaucus, 
            seatshare.change.countryFE.nocaucus, 
            seatshare.change.partyFE.nocaucus,
            type = "html",
            title = "Table B8. Negative binomial regression of post-election cabinet changes",
            out = "./tableB8.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Seat share change", "Total ministers", "Unavailable ministers", "Median tenure (current post)", 
                                 "Median tenure (any post)", "PM's first election", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(B25)", "(B26)", "(B27)", "(B28)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.change.bivariate.nocaucus), nobs(seatshare.change.controls.nocaucus), 
                               nobs(seatshare.change.countryFE.nocaucus), nobs(seatshare.change.partyFE.nocaucus)),
                             c("Pseudo R2", round(seatshare.change.bivariate.nocaucus.R2, digits = 3),
                               round(seatshare.change.controls.nocaucus.R2, digits = 3),
                               round(seatshare.change.countryFE.nocaucus.R2, digits = 3),
                               round(seatshare.change.partyFE.nocaucus.R2, digits = 3))),
            font.size = "small")
  
########################################################################################## 
  
  rm(list=ls())
  